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ABSTRACT 

The structure and behavior of barotropi cal ly unstable 
and stable wave disturbances in the vicinity of a zonally 
varying easterly jet are studied numerically with a 
linearized barotropic vorticity equation on a 3-plane. The 
easterly jet is approximated by a Bickley jet with a slow 
zonal variation. The numerical results are also compared 
with a simple me chan i Stic analytical model using the local 
phase speed and growth rate concepts. The results are 
grossly similar in several respects to that expected from 
the parallel flow theory of barotropic instability, however, 
the resultant structure of the waves causes a spatial growth 
rate greater than predicted by the local growth rates 
computed wfth a parallel flow model. In the stable region, 
the structure leads to strong dynamic damping. When a 
uniform advective velocity is added to a variable mean flow, 
the differences between the behavior of the computed waves 
and that implied by the parallel flow theory are somewhat 
reduced. The waves remove kinetic energy from the mean 
flow and most of this energy is removed on the downwind side 
of the jet. The computed structure and behavior of the waves 
have a number of features that resemble those observed in 
the vicinity of the upper troposphere easterly jet during 
the summer monsoon. 
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I. INTRODUCTION 



In the Northern Hemisphere summer, a moderately strong 
easterly jet exists south of the Tibetan high near the 200 mb 
level (Kri shnamurti , 1971a, 1971b). Synoptic-scale moving 
disturbances occur at the level of the jet, and it appears 
that these disturbances arise from barotropic instability of 
the mean flow. The jet contains regions of large vorticity 
gradients where the necessary condition for barotropic 
instability is sometimes locally satisfied. If the observed 
disturbances were the result of barotropic instability, they 
would extract energy from the mean zonal flow and the plane- 
tary-scale waves, since the latter combine with the zonal 
flow to give the large vorticity gradients south of the 
Tibetan high. In fact, Kanamitsu e_t al_. ( 1 972 ) have shown 
that wavenumbers 6-8 in the wind spectrum in the region 
between 15S and 15N receive energy through barotropic 
interaction with the zonal and wavenumber 1 flow. 

Furthermore, Krishnamurti (1971a) studied the wavenumber 
spectra of the meridional wind component for selected 
tropical latitudes and he found a peak in the spectra near 
wavenumbers 6-8 at latitudes near the easterly jet. 

Colton (1973) studied barotropic interactions between 
quasi -stati onary long waves and transient synoptic waves 
using a semi -spectral numerical model. His long waves were 
forced with a specified divergence field following the 
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diagnostic model of Holton and Colton (1972). Results from 
a long term integration of his model reproduced several 
features of the upper tropospheric general circulation. In 
particular, his model simulated the disturbances observed by 
Krishnamurti (1971a, 1971b) in the vicinity of the easterly 
jet. His model showed that wave disturbances entering 
upstream of the easterly jet regime grow and increase their 
speed and wavelength as they approach the longitude of the 
jet velocity maximum, achieve their maximum amplitude a 
considerable distance downstream of the jet maximum, and 
eventually decay as the disturbances move out of the easterly 
jet regime. Furthermore, the high and low pressure centers 
of these transient wave disturbances occur on the wings of 
the central axis of the jet stream. Colton concluded that 
the transient disturbances are due to scale interactions 
tnvolvtng short-term barotropic instability. 

In both Kri s hnamurti 1 s observational and Colton's 
numerical studies, the zonal variation of the jet apparently 
has significant effects on the dynamic behavior of the 
transient disturbances. This study is an attempt to under- 
stand better such effects. This is a linear study, therefore, 
the re suits can only be suggestive of the complicated nonlinear 
behavior of the real atmosphere. However, since the ampli- 
tudes of synoptic waves are generally of smaller magnitude 
in the tropics than in middle and higher latitudes, a linear 
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stability approach should be a reasonable approximation to 
the real atmosphere. 

This study is an extension of the linear barotropic 
instability theory developed over the years by various 
investigators. The extension of this study is the zonal 
variation of the basic flow. .Rayleigh (1890, 1913) developed 
the concept of hydrodynamic instability and Kuo (1949) 
extended this concept to a rotating atmosphere by including 
the beta term. Since then, many investigators have studied 
the stability of barotropic zonal flows. Kuo (1949, 1951), 
Lipps (1962), Yanai and Nitta (1968) investigated barotropic 
instability of parallel symmetric westerly jet flows. Mitta 
and Yanai (1969), Yamasaki and Wada (1972), and Kuo (1973) 
extended this theory to parallel symmetric easterly jet flows. 
Lipps (1965, 1970) and Kuo (1973) further extended the concept 
to parallel asymmetric barotropic zonal flows. In general, 
when the necessary criterion for linear barotropic instability 
is satisfied, normal mode solutions for barotropic unstable 
waves have certain properties: 1) they exist within a range of 
intermediate wavelengths, 2) they have a latitudinal tilt 
opposite to the shear of the barotropic zonal wind and this 
tilt creates momentum fluxes which redistribute the kinetic 
energy from the mean zonal flow to the disturbance, 3) their 
phase speeds in a westerly jet are less then the speed of the 
maximum westerly wind, and 4) their phase speed in an 



easterly jet may be greater than the maximum easterly wind 
(Pedlosky, 1964; Yamasaki and Wada, 1972). In this case, 
however, the growth rates are generally very small. 

Nitta and Yanai (1969) modified the concept that 
unstable solutions exist only within an intermediate range 
of wavelengths. For an easterly sinusoidal jet flow, they 
found a distinct short wave cut-off but no apparent long wave 
cut-off for instability. Yamasaki and Wada (1972) modified 
this by showing that the long wave cut-off is dependent on 
the strength and sharpness of this sinusoidal easterly jet. 
For very strong and sharp velocity profiles, the long wave 
cut-off approaches infinity. For relatively weak but still 
unstable profiles, a finite long wave cut-off exists for 
this jet. Kuo (1973) referred to this long wave cut-off 
region as the modified Rossby regime in his numerical study 
for an easterly Bickley jet. 

All these stability studies have used parallel baro- 
tropic flows, i.e., there is no zonal variation in the basic 
flow. Lorenz (1972), however, investigated primarily by 
analytical means the barotropic instability of a flow pattern 
which varies with longitude. The basic flow is a neutral 
Rossby wave superposed on a uniform westerly flow. Essen- 
tially, this flow depicts the progression of large scale 
waves embedded in a westerly current. Zonal flows of mid- 
latitudes are generally considered to be barotropi cal ly 
stable, but Lorenz showed that a uniform zonal flow together 
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with superposed neutral Rossby waves may be unstable with 
respect to further perturbations. Lorenz found that the 
growth rate of the perturbations is comparable to the growth 
rate of the errors of large numerical atmospheric models. 

Based on this, Lorenz suggested that barotropic instability 
may be partly responsible for the unpredictability of the 
real atmosphere. 

In addition to the regular or normal mode eigensolutions, 
there are "singular" or "continuum" mode solutions which 
have been discussed by Rayleigh (1913)., Case (1 960), Pedlosky 
(1964), and Yanai and Nitta (1968). These "singular" 
solutions correspond to continuous eigenvalues of phase 
velocity (c) which are equal to the basic flow, U'(y), some- 
where in the zonal current. Since these "singular" solutions 
are continuous, there is an infinite number of solutions. 

Case showed that these continuum modes are needed to form a 
complete set of solutions. He also showed that the disturb- 
ance formed by the sum of the continuum modes has a 
y-structure tilt in the same sense as the basic flow shear 
and that it usually decays as 1/t or faster, where t repre- 
sents time. 

It is clear that in the real atmosphere, "mean flows" 
vary both in space and time. They are neither purely 
barotropic nor purely baroclinic, and in including both 
effects in a linear stability study is a very difficult task. 
Most of the investigators have examined the dynamic 
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Instability problem by studying each effect separately. 

This study, however, focuses only on barotropic instability. 
We hope that a better understanding of the kind of upper 
tropospheric waves studied by Krishnamurti (1971a, 1971b) 
and Colton (1973) near the easterly jet can be achieved by 
studying the behavior of barotropic waves in a region of a 
variable mean easterly wind. We are especially interested 
in the situation where the waves move into and out of baro- 
tropically unstable regions. Therefore, the objective of 
this work is to determine the dynamic stability behavior, 
the structure and the energetics of barotropic waves propa- 
gating in a zonally varying mean wind. 

With a very simple analytical model, Pedlosky (1976) 
studied the finite amplitude dynamics in a zonally varying 
baroclinic current. The mean flow of his model changes 
abruptly from a weakly unstable regime with a constant growth 
rate to a stable regime downstream. He found that disturb- 
ances may propagate into stable regions substantially 
undiminished, retaining a considerable memory of its history 
in a locally unstable region. 

The mean flow in our barotropic study, however, has a 
full downstream variation such that the stability of the 
mean flow ranges from locally strong instability to locally 
strong stability. 
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Our approach is as follows: firstly, we develop a 
numerical model based on the linear non -di vergent barotropic 
vorticity equation on a beta-plane incorporating the dynamics 
of transient barotropic waves in a region of variable mean 
wind. The mean zonal wind is an easterly hyperbolic secant- 
squarred ( B i ck 1 ey ). jet and the mean meridional wind is derived 
in such a way that the mean flow is non-di vergent. Secondly, 
using this numerical model, we conduct selected case studies 
to determine fundamental dynamic stability and energetic 
properties of these barotropic waves. Thirdly, we compare 
the results obtained from the numerical model with the 
parallel flow theory. This is done by constructing a simple 
mechanistic analytical model which incorporates the local 
stability concept of the parallel flow theory. This compari- 
son gives us further insights on the effect of the variable 
mean wind. 

The domain of the numerical model consists of an open 
channel with rigid walls at the north and south boundaries. 
The basic flow is a slowly varying easterly Bickley jet. 

A periodic forcing is applied on the inflow or eastern 
boundary such that periodic perturbations are generated from 
this boundary into the channel flow which represents waves 
moving into the region from the east. The western boundary 
condition allows the waves to move out of the channel. As 
the waves move through the region, they grow or decay in 
relation to the local stability properties of the mean flow, 
whereas at each point the fields vary periodically. 
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II. NUMERICAL MODEL 



A. GENERAL FORMULATION 

The governing equations are the barotropic vorticity 
equation using a beta ( 6 ) plane approx imation: 

If + u H + v !f + ev ■ ■ v < 2 -’) 

and the non-di vergen t continuity equation 



3 U 3 V 
3 X 3y 



0 , 



( 2 . 2 ) 



where 



3 v 3 u 
3x " 3y 



(2.3) 



Here 8 is the north-south gradient of the earth's vorticity 
and is given by its value at 10° latitude. Q is a forcing 
function representing non -ba rotropi c effects which is 
required to maintain the vorticity field, is a frictional 
coefficient. Since the flow is two-dimensional and non- 
di vergent, a stream function ( 4 )) is de fined such that 



u = - 



3 ip 

ay 



_ 3 \jj 

3 X 



(2.4) 



thus 



? = 7 I 



(2.5) 



Here u and v are the velocities of the flow in the x and y 
direction, respectively. Equation (2.1) now becomes 

-L v 2 * - f± iiii ♦ e f± - q - o.v 2 * (2.6) 

3 1 v 3y 3X 3 x 3y 3X f 
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The governing equations are linearized by separating 
quantities into the following form: 

*(x»y»t) = $(x,y) + ^ ' ( x ,y , t ) . (2.7) 



Here the bars are the basic state variables and the primes 
are the perturbation quantities. The mean s treamf uncti on 
i)(x,y) satisfies the mean vorticity equation when a non- 
barotropic source term is added. The basic state equation 
is therefore 



ii + ii nil + oil - 5 

ay ax ax ay p ax y 



D f v 2 ijj 



(2.8) 



which is non-linear. 0 is needed to maintain the basic 

vorticity field in steady state. The resultant linearized 

perturbation equation is: 

- 2 - 2 1 2 - 
a 2 , _ a av \p 1 + a_^ 97 ^ 1 _ 341 av f 

at ^ " ay ax ax ay " ay ax 



+ 111 

ax 



nH + o ni 

ay p ax 



D f vV 



(2.9) 



The advective terms in (2.9) can be written in Jacobian form 
such that 



V 3 ^ - - J(4j,V%')- 0(^' ,V^^j) - 8 • - D^V^' , (2.10) 

where (2.10) is a Poisson equation with the tendency of the 
perturbation stream-function as the dependent variable. The 
right side of (2.10) is the forcing or non-homogeneous part 
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of the Poisson problem for 3ij>'/3t. The north-south boundaries 
are rigid walls where 



= 0 , • 



( 2 . 11 ) 



and 




( 2 . 12 ) 



The inflow or eastern boundary conditions are specified such 
that 



where to is the specified frequency of the forcing, A(y) 
and B Cy ) are the y-structure coefficients of the forcing 
which are described in Chapter IV. The purpose of this 
forcing is to introduce waves into the region from the east. 
The other condition on the inflow boundary is 

2 

C ' (0 ,y ,t) = 7%' = - kV + . (2.15) 

ay 

This expression follows because it is expected that i p' will 

i k x 

have a spatial variation of the form e . The determination 
of k is discussed later. This boundary forcing is expected 
to lead to periodic wave disturbances throughout the domain 
after a certain time period of numerical integration. It 
turns out that the outflow boundary conditions are crucial 



4»'(0,y,t) = A (y ) sin(wt) + s(y) cos(wt), 



(2.13) 



and 



||— ( 0 ,y , t ) = A(y )cocos (ut ) - 0 (y )us i n (ut ) , 



(2.14) 
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for attaining this periodic state, because the wave disturb- 
ances must be able to propagate freely out of this boundary. 
Therefore, a Sommerfeld radiation condition (Pearson, 1974) 
for both the tendency and the vorticity is used to approxi- 
mate this mechanism on the outflow boundary: 

ft <lf) - - ‘r fir (If). < 2 - 16 ) 



iil = . c Hi 
dt r ax 

o 



(2.17) 



Here c r is a specified constant phase velocity. Pearson 
o 

(1974) showed that the Sommerfeld radiation condition is a 
consistent boundary condition for numerical models of initial 
value systems admitting dispersive waves. The fundamental 
problem with most outflow boundary conditions is the reflec- 
tion of incident waves from a boundary back into the interior 
region. This is usually disastrous for numerical models 

admitting dispersive waves. Pearson showed that, if c is 

r o 

chosen judiciously, the longwaves move smoothly through the 
boundary. For short wavelengths, however, an area with a 
large coefficient of viscosity near the boundary is often 
required to help control the reflection problem. According 
to Pearson, the amount of damping is proportional to the 
wavenumber. A Matsuno finite difference scheme is used in 
our numerical model. This scheme has a tendency to damp the 
short waves thus we need only to use a relatively small 
coefficient of friction to control wave reflection at the 
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outflow boundary. A friction coefficient of D^. = .15x10 
sec”^ is used which is equivalent to an e-folding decay time 
of approximately 7.7 days. 

B. ZONALLY VARYING BASIC FLOW 

The basic velocity field is an easterly Bickley jet 
defined by 

u(x,y) = - U ( x ) sech 2 (j^y) - U Q = - ||- . (2.18) 

Here d(x) is a characteristic length scale of the jet and is 
related to the half width d(x) (Kuo, 1973) by 

d ( x ) = - 1.76 d(x) . (2.19) 



U(x) is the velocity of the jet at y=0 and U Q is a constant 
velocity. The basic flow s treamf un cti on (i) is specified 
to be constant at both the southern ( y = - D ) and northern (y = D) 
boundaries as follows: 



*(-D) = 

( D ) - C 2 



, where (2.20) 



y 

£(x,y) = u(x) J sech 2 ( dT^y) dy + u 0 f dy + ^ _I 

-D -D 



( 2.21 ) 
y 

i ^ ( y ^ Hw 4- II C Hu 4- i” ( _ Q ) } 



and 



D 



$(d) = u ( x ) J sech 2 (dTxT^ dy + u o I dy 



+ $(-D) = C 2 . (2.22) 
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For a symmetric jet profile about y = 0, C ] = C 2 . Equations 
(2.21 ) and (2.22) are conveniently simplified: 



<Hx,y) = U ( x ) d ( x ) {tanh(^-^y) + tanh(j^-y)} 



+ V + U o D + 



(2.23) 



and 



$(D) = 2U ( x ) d ( x ) tanh(j^J-) + 2U Q D + $(-D) . (2.24) 



From (2.23), the jet velocity at y=0 is given by: 



U(x) = 



$(D) - 2U 0 D - $(-D) 



coth (g^) . 



(2.25) 



2dTxT 



Therefore, if d(x) varies slowly in x, then, so does the 
basic flow. The x-variation for the characteristic length 
scale is given by 



Here x Q is the longitude where the x-variation of the cosine 
function starts and L is the wavelength of this variation. 

Figures 1-3 show the basic fields of streamfuncti on 
^(x,y), zonal velocity u(x,y) and vorticity £(x,y), 
respectively. Here L is set to 43,000 km. The domain is 
40,125 km long and 4,000 km wide (D = 2,000 km). At y=0, 
the basic zonal velocity is -15.51 m sec”^ at the inflow 
boundary (x = 14,625 km) increasing slowly to a maximum value 



d 



850 km + 350 km <cos[2ir 



1200 km 
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Fig. 1. Experiment I. The ^(x.y) field (x 10' nr sec 
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Experiment I. The u(x,y) field (x 10 m sec 
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Fig. 3 . Experiment I. The c(x,y) field (x 10 sec 



of -30 m sec"'' at x=0. From this longitude the central 
velocity slowly decreases downstream to a minimum value of 
-13.4 m sec"' at x = -21,375 km. Between this longitude and 
the outflow boundary (x = -25,500 km), the basic flow is 
parallel. From Eq . (2.26), we note that the characteristic 
length scale, d(x), varies between 500 km at u(x,o) v and 

CT1 a X 

1 200 km at u ( x , o) • . 

min 
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III_ LOCAL STABILITY OF BASIC FLOW 



The local growth rate of the variable mean flow of the 
numerical model is first determined by a parallel flow (i.e., 
no x-variation) numerical model of Williams e_t aj_. (1971) in 
order to gain some insight on the stability characteristics 
of the mean flow. This model is hereafter re fere need as the 
parallel flow model. By setting a ip / 3 x = 0 in (2.9) of the 
numerical model, the governing equation of this model is 



The symbols of (3.1) have the same meaning as the numerical 
model of Chapter II, except that the basic flow is given by 



where U is a specified constant that scales the magnitude of 
the central velocity of the Bickley jet (y=0). The character- 
istic length d is also a specified constant. Assuming that 
all perturbation quantities are periodic in x, Eq . (3.1) is 
finite Fourier transformed in x with wavenumber k, and is 
solved as an initial value problem. This approach gives the 
phase speed, growth rate, and the wave structure of the most 
unstable mode. It is convenient to write 4 > ' in the following 
form : 




(3.1) 



u(y) = - U sech 2 (<j) , 



(3.2) 



^'(x,y,t) = A(y , t ) cos(kx) + B(y,t) sin(kx). (3.3) 
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Equation (3.3) is substituted into (3.1), and the coefficients 
of coskx and sinkx are separately set equal to zero. The 
cos kx term gi ves : 



[ 




k 2 ]|A 
J at 



2 g- 

k C u - k 2 } + (6 - - — *■} ] B 

3y ay 

2 • (3.4) 

Df [~ — t “ ^ 2 ]A , 

3y 



and the sinkx term gives: 




k 2 l— 

J at 




- M 



ay 



k 2 } + (B 
- k 2 ] B . 



^|}]A 

ay 



(3.5) 



The boundary conditions are 



1 ( x , D , t ) = 1 (x,-D,t) = 0. 



(3.6) 



Equations (3.4) and (3.5) are written in finite difference 
form such that the second derivative with respect to y of a 
typical variable A is approximated as 



cty - c ^ 1 

ay ,• 



- 2A .i + A .i-i 
Ay 2 



] , 



(3.7) 



where j is the y-grid index and Ay is the distance between 
grid points. The model has 32 grid intervals and Ay is set 
equal to 125 km (the width of the channel is 4000 km). This 
y-grid structure is the same one used in the complete 
numerical model. Centered time differences are used for all 
quantities except those involving friction. The time step 
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(At) is set to 1 hour. The friction terms are evaluated at 
the previous time step in order to insure linear computa- 
tional stability (Haltiner, 1971). The integration begins 
with a forward time step. The boundary conditions (Eq. 3.6) 
become 

A = B = 0 aty = +D . (3.8) 

Equations (3.4) and (3.5) are solved for the tendencies by 
the exact method of Richtmyer (1967). These equations can be 
solved numerically as a function of 8, D f , u ( y ) , and k for 
any initial conditions. These equations are integrated until 
their solution becomes exponential in behavior. The basic 
flow u(y) is specified for a number of selected longitudes 
of the variable mean flow of Fig. 2. The wavenumber k is 
also specified for each integration. 

In general, these equations have a set of discrete normal 
mode solutions as well as a continuous spectrum of solutions 
(Case, 1960, Pedlosky, 1964 and Yanai and Nitta, 1968). Only 
the normal mode solutions can give significant growth and the 
most unstable mode will dominate after a sufficient period 
of time. 

This parallel flow model is numerically integrated to 
150 days. The initial disturbance amplitude has a north- 
south structure with no tilt. The eigensolution obtained 
includes wave structure of the most unstable mode, and its 
growth rate and phase speed. By selectively testing different 
wavelengths, the most unstable wavelength is determined. 
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Figures 4-7 illustrate the eigensolutions of the most 
unstable discrete mode for the selected jet profiles. Since 
the jet flow and the wave disturbance are both symmetric 
about y= 0, only the lower half of the y domain is shown. 

It is interesting to note that the eigensolution amplitude 
fias three maxima, one at y = 0 and one on each wing (only lower 
wfng is shown) of the jet approximately 600-800 km from y=0. 
It is apparent that when the jet is relatively sharp (i.e., 
d < 700 km), the central maximum predominates. When the jet 
is relatively smooth and broad (i.e., d > 800 km), the 
maximum on the wings predominates. Figure 8 shows the 
growth rate (n) corresponding to the most unstable wave- 
lengths (L) as a function of x based on the parallel flow 
model. In a local stability sense, the most unstable 
wavelengths range between 3650 km at x=0 where the jet has 
its maximum velocity and 4600 km near the inflow and outflow 
regions of the jet regime. We observe that locally the 
greatest instability based on both tilt and growth rate is 
indicated where the jet achieves its maximum central velocity 
(u = -30 m sec" 1 ). This parallel flow model can only give 
the most unstable discrete mode; it can not depict any 
dynamic damping. This is reflected in Fig. 7(b) which 
corresponds to x = -22,500 km, outside the unstable region. 
Here the solution exhibits essentially no tilt and the growth 
rate asymptotically approaches the linear frictional damping 
rate. This is further discussed in Chapters VI and VII. 
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Fig. 4. Parallel flow model: wave structure of the most unstable mode 
for longitudes (a) x = 0, and x = +^3 750 km. 
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Fig. 5. Parallel flow model wave structure of the most unstable mode 
for longitudes: (a) x = +7500 km, and (b) x = +^9 375 km. 
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Fig. 6. Parallel flow model wave structure of the most unstable mode 
for longitudes (a) x = +11250 km, and (b) x = +13125 km. 
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Fig. 7. Parallel flow model wave structure of the most unstable mode 
for longitudes (a) x = +14625 km, and (b) x = -22500 km. 
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Fig. 8. (a) The most unstable wavelength L(x), and (b) corresponding 

local growth rate n (x), based on the parallel flow model. 



These results are consistent with those calculated by 
Kuo (1973) in a numerical study of a parallel Bickley jet 
and, in a local stability sense, may be viewed as a first 
approximation of the behavior of moving waves within a 
zonally-varying mean wind. 
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IV . FINITE DIFFERENCE EQUATIONS 



The governing Eq. (2.9) of the complete numerical model 
is written with the advective terms in Jacobian form: 

v 2 ||^- = - J($i,vV) - J(^',v 2 {p) - eff^- - D f vV , (4.1) 

where 



j(a ,b) = - iiii 

' ' ax ay ay ax 



(4.2) 



for any two scalar quantities a and b. The Arakawa (1966) 
finite difference approximation for the Jacobian is used. 

This scheme conserves both the mean square vorticity and the 
mean kinetic energy. 

We define the tendency of the disturbance s treamf un cti on 
by 



T? i 5 [If-] 1 , (4.3) 

1,J at 

where superscript t is the current time, and subscripts i 
and j refer to the x and y grid points, respectively. In the 
following finite difference equations. At is the time step 
and is set equal to one hour, ax and Ay are the x and y 
grid point intervals, respectively. Equation (4.1) is written 
using the Matsuno or Euler-backward finite difference scheme: 






- 6 



-1 ot-At 

. 1 >J 2ax D fV »{ tj ’ 



(4.4) 
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(4.5) 



* tt 

ii > : . = ii) . . + t . .At , 

^1 > J M , J 1 , J 



z * * * 

V T i>j = -X jCMV] -X - c ^ 1 , v 2 ^ ] 

i > j i > j 



- 6 






i> • 



1 + 1, j ' y i-K.i 

2 Ax 



- D 



__2 



t-At 



f V*i,j 



(4.6) 



, t + At _ , t T * . f 

*i,i ' * i .. :i + T i.j At 



(4.7) 



The finite difference form for the Laplacian of the tendency 
is given by: (4.8) 



yt _ o jt +T^ T~ ^ _ ?T^ + T** 

tp- T v _ 1+1, j i’ , .1 i , J + 1 . i,j + l d i » J i , J - 1 

V T i i = 2 2 

1,J AX Ay 



2,t 



and for the Jacobian terms, it is given by: 



© t _ , 

■JH .• " 4ixAy 
'* > J 



( ^i + l ,j " Vi (c i\.j + 1 ‘ ? i J-l } + 



^ ,j + l “ *i ,j-l } (c i+l ,j " ? i 



l t "] 

i - 1 , J ) J , 



(4.9) 



© t , 

JT. . " 4AXAy 

1 5 J 



[ ^ i + 1 , j ^ i + 1 , j + 1 " ? i + 1 , j - 1 ^ + 



" Vl ,j + l ' C i -1 , J -1 } + 



(4.10) 



, t 



. t 



' *1 ,j + l (? i + l ,j + 1 " c i-l ,j + 1 ) + 



+ *i ,j-1 (c i + l J-l " c i-1 J " 1 } 



■ t 4 
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and 



(D t 

JJ\ .. 



4AXAy 



i , J 






1 J + 1 v ^ i + 1 J + 1 ' v i-l ,j + l ) + 



- c 



i J-l ( *i + l J-l ' *1-1 J-l } 



(4.11) 



? i + 1 , j ^ * 1 + 1 » j + 1 " * 1 + 1 , j - 1 ^ + 



+ ©1 J ( *1-l J + 1 " *1-1 J-l 



where 



©t ©t 

JC = J . . + 5F. . + "3T 



t 



i » j i ,j 



J_JL 



(4.12) 



Jr . [t|/ ',v 3 is expanded in the same manner as for 

i J 



tt.’ 2 *'] • 
i J 



The Poisson Eq. (4.1) is solved for the tendency with 
a direct method developed by Sweet (1971). This direct method 
solves a finite difference approximation to Poisson's equation 
on a rectangular domain with Dirichlet boundary conditions. 

The finite difference form of the boundary conditions are as 
f o 1 1 ows : 

a) Northern (j=0) and southern (j=0) boundaries: 



ij©. = 0 : i=0,I for j = 0 and j=J, 

1 > J 



(4.13) 



4 



,t 

1 »o 




1 = 0 , 1 , 



(4.14) 
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(4.15) 



. t 

i ,0 



, t 

i , J - 1 ’ 



i = 0 , 1 



b) Inflow or eastern boundary (i=I, j=0, j) 



. = A. sin(uilAt) + B. cos(aj£At), 
^ J J J J 



Tr • = A .COCOS (o)2,At ) - B. ojsi n(cu£At) , 
1 J J J J 



and 



-k? . *:*. + -b-iili 

1,3 I»J 



. t 



1,3 






Ay 



(4.16) 

(4.17) 



(4.18) 



where to is the specified frequency and k T . is a specified 
wavenumber structure. Time is discretized by: 

t = £At, where i - 0,1, 2, 3, 4 ... (4.19) 

A,, and B. are the y-structure coefficients of the forcing 
j , J 

(4.14). These coefficients were determined from the eigen- 

solution of the parallel flow model. The velocity profile 

for determining A- and B. of the inflow boundary condition 

3 3 

was selected from the inflow region of the numerical model. 
This profile is described in Chapter V. 

c) Outflow or western boundary (i=0, j=0, J) 



j t-At T t 

ft - Q + Lj r 

0,j l+c At/AX 1+AX/(C At) 
r o r o 



(4.20) 



s 



, t 
0,3 



' t-At 

C 0,.i 

1+C At/AX 



i i 4 ... _ 

1 +AX/ ( C At) 



(4.21) 
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(4.22) 






. t 

0 , j 



. i t-A t T * 

' *0.j + T oj 



A t . 



The tendency . (4.20) and- (4.21) were derived from 

the analytical Sommerfeld radiation condition (2.16) and 
(2.17), respectively. Since (2.16) is not in the usual form 
of the Dirichlet boundary condition, the direct Poisson 
equation solver, POISDD, had to be modified. Details of this 
modification are given in the Appendix. 
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V. EXPERIMENT I 



The investigation of the numerical model of the variable 
mean wind comprises two experiments. Experiment I is the 
principal one and will be discussed in detail in this chapter. 
In Experiment II a constant velocity is added to the mean 
flow and its results will be discussed in Chapter VII. 

The forecast equation is integrated in time from an 
initial state of \p 1 = 0. Tne periodic forcing on the eastern 
boundary causes the interior streamfuncti on to grow and the 
integration is continued until the time variation is periodic 
everywhere with the forcing frequency. By day 70, the model 
has achieved the fully periodic state and the wave packet 
envelope has become q uas i -s tati onary . The forcing frequency 
is varied until the value which gives a maximum perturbation 
amplitude is found. 

The following values are used in Experiment I: 

2D = 4000 km, x^-x^ = 401 25 km, u ( 0 , 0 ) = - 30 m sec"^ , 

U Q = 0, L = 43000 km. Ax = 375 km, Ay = 125 km. 

For these parameters the maximum response occurs for a 
forcing period of 3.25 days. It requires several complete 
integrations to refine the inflow and outflow boundary 
conditions for the specified forcing period. Wavenumber (k) 
for the inflow vorticity boundary condition (4.13) was 
determined by observing the predominant wavelengths near the 
inflow region. The radiation phase velocity (c ) in the. 
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outflow boundary condition (4.15) was similarly determined 
by observing the predominant wavelengths near the outflow 
region. Since the domain is periodic, the phase velocity is 
approximated by the simple relation: 




Here u is the local frequency which, in the fully periodic 

state, is equal to the forcing frequency, and k Q is the wave- 

number measured near the outflow region. Thus, for Experiment 

I, a phase velocity (c ) of -13 m sec"^ is used for the 

r o 

outflow Sommerfeld radiation boundary condition. 

A. RESULTS 

The solution for the Case I experiment becomes fully 
periodic after 70 days. The > ' - f i e 1 d at t=70 days is shown 
in Fig. 9. An entire train of barotropic waves actually 
exists throughout the length of the channel, but the waves 
upstream of x = -4,500 km and in the outflow region are not 
shown in Fig. 9 because of their relatively small amplitude. 
The maximum wave amplitude occurs at x = -12,750 km and is 
2 and 4 orders of magnitude larger than at the jet maximum 
(x=0) and the inflow boundary (x = 14,625 km), respectively. 
Figure 10 shows the lower half of phase angle tilt 9*(x,y) 
of the wave disturbance y-structure for various values of 
longitude x. We note that the waves upstream of approximately 
x = -20,000 km tilt opposite to the mean wind shear, which is 
necessary for barotropic instability. In fact the tilt near 
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Fig. 9. Experiment I. The i|» ' field at t = 70 days. 
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Fig. 10. Experiment I. The phase angle tilt 8*(x,y) for longitudes (a) x = 13125 km 
(b) x = 5625 km, (c) x = 0, (d) x = - 5625 km, (e) x = -13125 km, (f) 
x = -19500 km, and (g) x = -24375 km. 



tlie inflow boundary shows relatively weak instability and 
this tilt slowly increases, reaching a maximum at x=0 where 
the jet velocity is maximum. Further downstream, the tilt 
slowly decreases and eventually reverses near x = -20,000 1cm. 
This behavior is consistent with the growth rates for the 
most unstable wavelengths shown in Fig. 8 which were computed 
using the parallel flow model. Near the outflow boundary, 
however, the tilt is reversed. This indicates dynamic 
stability or a flow of energy from the disturbance back to 
the mean flow. Therefore, we observe that the tilt of the 
wave disturbance qualitatively adjusts to the local stability 
of the mean flow (see Figs. 4-7). We recall that the parallel 
flow model can only solve for the most unstable discrete mode, 
thus, dynamic damping is not indicated in Fig. 7(b) for the 
outflow region. 

Figure 11 shows the envelope of the wave packet, 
evaluated at y = -750 km, where the disturbance amplitude is 
large. This envelope <i>'(x)> is obtained by recording the 
maximum and minimum tp 1 values that occur at each longitude 
over a period of ten days after the solution becomes fully 
periodic. If a larger time interval were chosen, the envelope 
would not change. Note that the maximum amplitude occurs in 
the area where the local growth rate becomes zero (see Fig. 8). 
The smoothness of the fields in Figs. 9 and 11 indicate that 
the radiation outflow boundary condition is working properly. 
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Experiment I. Wave packet envelope, <^'(x)>for 
the complete numerical model. 



The tp ' fields for days 83 through 86 are shown in Figs. 

12 through 15, respectively, to illustrate a time sequence 
of the ip 1 barotropic wave train. A periodicity of 3.25 days 
can be determined by examining the time series at each grid 
point (not shown). We also note from this sequence that there 
are three maxima of ip 1 in the wave structure. One maximum is 
at y=0 and another is on each wing of the jet near y = +750 km. 
The maximum on the wing clearly predominates. Figure 16 shows 
the variation of wavelength, L(x), for y=0 and y = +750 km 
corresponding to the latitudes of the observed three maxima 
of ill'. The disturbance wavelength is about 4100 km initially 
near the inflow boundary. The wavelength near y = +750 km 
is larger upstream and smaller downstream of x=0 than at y=0. 

At y = +750 km, the maximum wavelength of 5060 km occurs 
approximately 950 km upstream of x=0 while the minimum wave- 
length of 3600 km occurs near the outflow boundary. At y=0, 
the maximum wavelength of 4900 km occurs about 950 km down- 
stream of x=0 while the minimum wavelength of about 4000 km 
occurs near x = -13,500 km with the wavelength increasing 
slowly further downstream to the outflow boundary. The 
range of the wavelength is between approximately 3600 km and 
5000 km, which is nearly the same range of values obtained 
for the most unstable wavelengths using the parallel flow 
model in Chapter III. Figure 17 shows the disturbance phase 
speeds, c f *( x )> for latitudes y=0 and y = +750 km where 

c *(x) is obtained from 
r 
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Fig. 12. Experiment I. The ip ' field at t = 83 days. 
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Fig. 13. Experiment I. The ij> ' field at t = 84 days. 
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Fig. 14. Experiment I. The ' field at t = 85 days. 
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Fig. 15. Experiment I. The ip ' field at t = 86 days. 
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Experiment I. Wavelengths L(x), from the complete numerical 
model for latitudes (a) y = 0, and (b) y = -750 km. 
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Fig. 17. Experiment I. Phase velocities c r * ( x ) , from the complete 

numerical model for latitudes (a) y = 0, and (b) y = -750 km. 



(5.2) 



c r *(x) - ^ 



Since the entire domain is periodic with the forcing 
frequency, c r *(x) has the same basic behavior as L(x) 

[Fig. 16] for the correspondi ng latitudes. Namely, upstream 
of x=0 the disturbance phase velocity is faster on the wings 
of the jet than at y=0, and downstream of x=0 it is slower. 
This behavior of both L(x) and c r *(x) illustrates the tilt 
behavior of the waves as they move into and out of the baro- 
tropically unstable region of the variable mean wind. More 
results are presented in the next chapter. 

The parallel flow model is also employed to determine 
a reasonably good forcing function that would generate 
periodic wave disturbances from the eastern boundary. This 
is accomplished as follows: 

Equation (2.14) defines the analytical form of the 
periodic forcing which is applied as an inflow boundary 
condition in the complete numerical model, i.e., 



Equation (4.17) describes the finite difference approximation 
to this forcing. The coefficients A ( y ) and B(y) of (2.14) 
are obtained from the ei gensol uti on determined using the 
parallel flow model for the following Bickley jet profile: 



H- ( 0 ,y , t ) = A(y )wcos (cut) - B (y ) wsi n (u>t ) . 



u(y=0) 

d 

L 
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962 km . • 
4,600 km 



-1 
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This profile is representative of the mean flow near the 
inflow region. Figure 18 shows the y-structure of these 
coef f i ci ents . 

Two other forcing functions, which are independently 
constructed from the parallel flow model, are also tested 
These functi ons are : 



fl” (0,y,t) = - u cos ( — — q^-) si n (ut ) , (5.3) 

and 

g 

g (0,y ,t) = - to { cos ( } sin(a)t), (5.4) 

where D is the half width of the channel (0 = 2000 km). 

The behavior of the wave disturbances for all three 
forcing functions are essentially the same except in the 
vicinity of the inflow boundary. The parallel flow solutions 
for the most unstable wavelength show a weak decay rate in 
the inflow region (see Fig. 8). This is because the linear 
friction decay rate is slightly larger than the weak local 
growth rate. We found that this weak decay behavior near 
the inflow boundary is achieved in the numerical model using 
forcing function (2.14). Use of the other two forcing 
functions, (5.3) and (5.4), realizes a weak secondary <\p '> m , v 
near the inflow boundary. Therefore, we found through use 
of the parallel flow model that the first forcing function 
provided the better inflow boundary condition than the other 
two, although no substantial difference occurs even if the 
other forcing functions are used. 
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Fig. 18. The coe-Pfic.ie.rrts A (y) and B.(y) 
of the eastern boundary forcing 
function used in the complete 
nume ri cal model. 
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B. ENERGETICS 



To understand the dynamics of barotropic instability, 
it is helpful to derive the disturbance kinetic energy 
equation. A finite difference approximation form of this 
equation is used to compute the energy budget of the waves. 

The linearized momentum equations are: 



i “1 + + u'2i + vl^l + v'2^ + 



3 t 



3 X 



vii' 
3 X ' 3y 



ay 



- syv ' = - r If 1 - D f u '- 



and 



+ Qivl + u'li + v^- + v'|i + 
3 1 3X 3 X 3y 3y 



Syu 



.1 . 



- ay 

p 



D r v ' 



(5.5) 



(5.6) 



After multiplying (5.5) by u' and (5.6) by v' and adding the 
resultant equations, we obtain the disturbance kinetic energy 
equation: 



i*I + + v— 

3 1 3 X 3y 



= - u 



' 2 ii 

3 X 



. - 3 U i i 3 V 

-u'v 1 — - u'v' — + 



ay 



2 - 



(5.7) 



- v 



3 V 1_ 3 ( U 1 p 1 ) _ 1_ 3 ( V 1 P 1 ) 



ay 



3 x 



ay 



2D ok 



where 



,2 



k 1 s 



v ' 



(5.8) 



The disturbance fields are considered to have a time 
dependence of the form sin[wt + 0(x,y)]. Thus, when (5.7) 
is averaged over one period and over one latitudinal domain 
(-D<y<_D), the following energy balance equation is obtained: 
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(5.9) 



6 



7 



9 





0 . 



Term 1 is the advection of k' by the mean flow (u). Terms 
2, 3, 4 and 5 are the barotropic exchange terms with the 
mean flow. Term 6 represents the work done by the pressure 
field. Term 7 is the dissipation term and represents an 
energy sink. With the exception of term 6, all the terms i r\ 
(5.9) are determinable from the $ and $ ' fields. Thus, the 
pressure work term is determined as a residual in Eq . (5.9). 
The energy equation is evaluated after the complete numerical 
model achieves the fully periodic state. 

Figure 19 shows the energy balance as function of 
longitude. Clearly, the maximum value of each of the terms 
in the energy equation occurs downstream of x=0. Curve (a) 
represents the Reynold's stress term, <( -u 1 v 1 3 u/ 3y ) > , where 
the maximum value occurs approximately at x = -10,975 km. 
Curve (b) represents the advection term, <(-u3k'/3x)> , of 
k 1 by the zonal current. Curve (c) represents the linear 
dissipation term <(2D f k')>. We observe that the largest 
dissipation takes place near the maximum amplitude of the 
wave packet envelope, <$ 1 ( x ) > . Curve (d) represents the 

I It U A 



pressure work 




It also has a maximum 



near <$ 1 ( x )> max • The remaining barotropic terms are not 
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Experiment I. Disturbance energy balance as a function of longitude 
(a) Reynold's stress term <( u ' v ' 3 u/ 3y )>; (b) energy transport term, 
COuak ' /3x )>; (c) dissipation term <( 2D^k ' )>; and (d) pressure work 

term <(- ~ ( u * p ' )>. 



included in Fig. 18 because they are 2 to 3 orders of 

magnitude smaller than <(- u'v'3u/3y)>. 

Recall that for Experiment I, <a>'(x)> occurs approxi- 

max rr 

mately at x = -12,750 km. We observe that the advection of 

k' by u plays a dominant role in the energy balance. Once 

the fully periodic state is achieved, the generation of 

disturbance energy by barotropic instability essentially 

occurs in the region between the <\p'> (x = -1 2750 km) and 

max 

the jet maximum (x=0). It is in this region that the 
strongest positive correlations occur between u'v' and 
- 3 u / 3 y in the Reynold's term. This is consistent with the 
observed behavior of the waves (see Figs. 10 and 11) since 
this strong correlation of u'v 1 is dependent on maximum 
amplitude of ^'squarred and the maximum tilt of the wave 
structure with the shear of the mean flow -3u/3y. It is 
also from this source re g ion that u transports a large part 
of this generated k' downstream to a point centered near 
x = -15,000 km . 

Strongest dissipation of k' occurs in the region centered 

near <^'(x)> . It is in this region that the strongest 

m ci x 

spatial gradients of ip ' occur. The largest values of pressure 
work also occur in this region where the pressure work term 
is computed as a residual. We note that its location of 
occurrence is consistent with the large amplitudes and spatial 
gradients of ip 1 of the region. We also note that <(-u3k'/3x)> 
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vanishes at <tJ» ' (*.) > max since k' is maximum at this longitude. 
The disturbance energy balance which is shown in Fig. 19 
graphically illustrates the basic energetics involved in 
supporting the large amplitudes and spatial gradients of ip ' 
a considerable distance downstream of the jet maximum. 

Clearly, the Reynold's stress term, < ( u ' v 1 a.u / a y ) > , is the only 
dominant source for disturbance kinetic energy. 
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VI. SIMPLE MECHANISTIC ANALYTICAL MODEL 



A. FORMULATION 

In this chapter we develop a simple mechanistic analyti- 
cal model which uses the locally determined parallel flow 
solutions. The results of this simple model can be compared 
wfth those of the complete numerical model to obtain some 
measure of the adjustment of the transient barotropic waves 
to the local stability of the variable mean wind. 

The following equation allows for propagation and growth 
or decay: 



where rp‘ represents the disturbance streamfunction. c^ (x ) 
is the local phase velocity and n(x) is the local .growth rate. 
We observe that if c r and n are independent of x, then 
Eq. (6.1) is exact, and, in general, should give a reasonable 
approximation to the downstream variation of tp ' . 

Let us consider the case where ip ' has a solution of the 

form : 



where oj i s a specified frequency. After substituting (6.2) 
into (6.1), we obtain 



It- + c r (x) = n(x) * 



( 6 . 1 ) 




( 6 . 2 ) 



iwF' + c r (x) = n (x ) F ' . 



(6.3) 
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Rearranging (6.3), we get 



1 dF 1 



- 1 to 



(6.4) 



F 1 dx 



c r (x) 



which has a solution of the form of 

x 




(6.5) 




( 6 . 6 ) 



F* (x) = Fq exp 



Some fundamental insight can be gleaned from (6.6). Here 
the amplitude of 1 must be specified at some initial point, 
x Q , which could be the inflow point. The local wavenumber 
is to/ Cp and the spatial growth rate is n/c r . These are the 
real and imaginary parts of the wavenumber, if we were to 
write 1 = A exp (ikx-iwt). 

In this simple analytical model, we are interested in 
the growth and decay of the wave packet. If the oscillatory 
behavior of Eq. (6.6) is dropped, the following equation 
expresses the envelope of the wave packet: 



Equation (6.7) gives the exponential growth/decay behavior 
of a wave disturbance as it travels within its wave packet 
envelope. This behavior is determined by the integral 
effects of the local stability properties of the mean flow. 



x 




0 



(6.7 
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B. APPLICATIONS 

Figure 20 compares the envelope for tf/'(x) from the 
complete numerical solution (Fig. 11) with the F(x) given 
by Eq. (6.7) for two values of x Q . In this figure the 
amplitude is plotted against upstream di stance from the 
initial point x . The analytical model uses the locally 
determined parallel flow growth rates (n) and phase speeds 
(c r ) based on the wavelengths measured from the complete 
numerical solutions. Therefore, Eq. (6.7) is applied at 
latitude y = +. 750 km where the wavelengths are measured. 
The curve for F(x) can be shifted up or down by changing 
x q , but its shape does not change. In Figure 20 all the 
curves have a maximum amplitude at x= - 12,750 km where 
the local growth rate is zero; however, the simple model 
has a lower maximum growth rate and slower damping rates 
than the complete numerical model. 

The lower portion of Fig. 21 contains the wavelength 
measured from the numerical solution at y = - 750 km. In 
the upper portion of the diagram are the phase velocities, 
c from the parallel flow model and c^* from the complete 
numerical model for the same latitude. Here c r is computed 

l 

with the use of the wavel ength; shown . The two phase speeds 
have similar behavior although c^* is shifted slightly 
upstream and has larger variations. We also note that c^* 
is generally smaller than c r except for the region 
-1500 km < x < 6750 km. 
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Fig. 20. Experiment I. The wave packet envelopes. F (x) is from the 
analytical model for initial point (a) x 0 = -375 km, F|j ( x ) 
is from the analytical model for initial point (b) x 0 ? 1875 km, 
and <•/>*( x )> is from the complete numerical model. 
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Similarly, Fig. 22 contains the wavelengths and phase 

speeds for latitude y=0. Again, the two phase speeds have 

similar behavior, but at this latitude c * is shifted 

r 

slightly downstream and the differences between phase speeds 
is generally less than at y = -750 km. Again we note that 
c r * is generally smaller than c r except for the region: 

-2250 km < x < 2250 km. This indicates that both the wave- 
length L(x) and the phase velocity c r *(x) between latitudes 
y=0 and y = -750 km are maximum where the mean flow is 
strongest . 

Figure 23 contains the spatial growth rates, m and m*, 
from the parallel flow model and from the complete numerical 
model, respectively. For the parallel flow model, we obtain 

m = - n/c r . (6.8) 

The value of m* is computed directly from the envelope. We 
note in Fig. 23 that m* has a larger maximum than m and the 
maximum is shifted slightly downstream from the jet maximum. 
Both curves pass through 0 at x = -12,750 km which is the 
maximum for both wave packet envelopes. However, m* 
shows a much larger damping in the outflow region which can 
also be seen in Fig. 20. In fact the parallel flow solution 
damps at the rate given by the frictional coefficient 
divided by the phase speed. This is expected in the baro- 
tropi cal ly stable outflow region when using the parallel 
flow model because this model can only give the eigen- 
solution of the most unstable discrete mode. 



70 



LATITUDE: 



<- (W>l 000 L ) 1 



O CO ^ C\J O coco 

• •• •• 
uo coco 




i i ii iii 

( ^_33S Ui) 



71 



-18 -15 -12 - 9 - 6-303 691^ 

X ( 1 0Q0 KM) + 

Fig. 22. Experiment I. Latitude y = 0. (a) phase velocities, 

c r (x) from parallel flow model, and c r *(x) from complete 
numerical model. (b) Wavelength L(x). 
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-18 -1 5 -12 -9 -6 -3 0 3 6 9 12 

X ( 1 000 KM) + 

Fig. 23. Experiment I. Spatial growth rates, m(x) from the parallel 
flow model, and m*(x) from the complete numerical model. 



The solutions observed in the numerical model (see 
Fig. 10) toward the outflow boundary tilt in the same sense 
as the mean wind shear which gives dynamic damping. This 
tilt behavior is also indicated by the divergence of phase 
speed in the outflow region between y=0 and y = -750 km 
(see Figs. 21 and 22). This dynamic damping appears to 
be due to continuous spectrum solutions (Case, 1960; Yanai 
and Nitta, 1968) which are not included in the normal mode 
solutions that are employed in the simple integral 
[Eq. (6.6)]. In fact we observe that, in the outflow region, 
the dynamic damping is 1/t or faster, increasing near the 
outflow boundary. In Fig. 23 the m* curve is skewed slightly 
to the left with respect to m and the jet maximum. This 
could be the result of the tilt structure in the wave lagging 
spatially behind the expected value from the local stability 
conditions. This would give a smaller growth rate on the 
upwind side of the jet maximum and a larger growth rate on 
the downwind side. This effect is indicated by the asymmetry 
of the m* curve with respect to the m curve in the unstable 
region of the mean wind. However, this is not the only effect 
involved because the most striking feature of Fig. 23 is 
the fact that the maximum value of m* is significantly 
larger than the maximum value of m. We have already shown 
that the most important disturbance energy production term 
(see Fig. 18) in the energy computations with the numerical 
solutions is proportional to 
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< ( -u ' v 1 9u/ 9y ) > . 



This term is the only source term for the parallel flow 
model, and it depends on the phase tilt in the disturbance 
field. Figure 24 compares the phase tilt of the numerical 
model solution, 9 *, and that of the most unstable wavelength 
of the parallel flow model, 9 , at x = 0, x = 3750 km, and 
x = 13125 km. It is evident that at the jet maximum the 
tilt from the numerical model is significantly larger than 
for the tilt from the parallel flow model. This is also 
evident but to a lesser degree at x = 3750 km. These results 
are consistent with the larger growth rates for the variable 
jet flow. This stronger tilt is also indicated by the 
comparison of the c^* and c r curves in the region near and 
slightly upstream of the jet maximum for latitudes y=0 and 
y = -750 km (see Figs. 21 and 22). We note that at 
y = -750 km, c r * is relatively larger than c r near the jet 
maximum. This difference is significantly reduced and 
shifted slightly downstream at y=0. We also note that near 
the inflow region (x = 13125 km), the phase tilt of the 
numerical solution lags slightly behind the parallel flow 
solution (see Fig. 24c). Apparently the downstream varia- 
tion of u augments the phase tilt which gives a larger 
growth rate. 

In comparing the parallel flow theory and the complete 
numerical model, a source of error could be introduced 
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Fig. 24. Experiment I. Comparison of phase tilt of the complete 
numerical model, e*, and of the most unstable wavelength 
of the parallel flow model, 0 q, for longitudes (a) x = 0 
(b) x = 3750 km, and (c) x = 13125 km. 



through the zonal resolution. This is because in the former 
the zonal variation is specified by a spectral representa- 
tion, while in the latter it is resolved by finite differences. 
Since the wavelength appearing in the complete numerical model 
varies between 10 and 13.5 ax, such resolution normally gives 
a good approximation to the exact solution (Haltiner, 1971). 
Thus, this error should not be very significant. In order to 
confirm this we carried out an experiment in which a parallel 
mean flow is specified by the full two-dimensional model 
ftntte differences of the complete numerical model. The 
growth rate and phase speed calculated from the experiment 
agree very well with those obtained by the semi -spectra 1 
parallel flow model, the difference being only 1-2%. Thus 
we may conclude that the differences between the numerical 
solution and the local solution of the parallel flow theory 
are genuine and are not due to the differences in the model 
resol uti on . 

C. EFFECTS OF LINEAR FRICTION 

Figure 23 shows that the local growth rate, determined 
by the parallel flow model, changes sign at x = -12,750 km 
and that it approaches the linear friction decay rate near 
the outflow region. We also note that the maximum aplitude 
of both the analytical and numerical model wave packet 
envelopes occurs at the lontitude where n(x) vanishes. 

Based on these results, <4' 1 > Iliax is expected not only to 
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decrease in magnitude but to occur further upstream if the 
linear frictional rate is increased. Likewise, the opposite 
effect would be expected if this friction were decreased. 

This hypothesis is tested by rerunning Experiment I using 

-5 -1 

a larger friction of = .25x10 sec which is equiva- 
lent to e-folding decay time of 4.63 days. We found that 
<< ,,> max“ occurred a PP roxlmate ^ 1 700 km further upstream and 
with an amplitude of an order of magnitude smaller (.4 vice 
5.9) than for the original experiment. 
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VII. EXPERIMENT II 



In Experiment II a constant advective speed of 
Uo = -5 m sec"^ is added to the mean flow, otherwise the 
parameters are the same as those used in Experiment I. We 
hope this modification will shed some light on the adjustment 
process of a faster moving wave chain. As in Experiment I, 
the analytical model also uses locally determined parallel 
flow growth rate (n) and the phase speed (c) based on the 
wavelength measured from the complete numerical solution. 

The local phase speed c is now related to the Doppler- 
shifted phase speed c by c = c + U o . The maximum \p 1 response 
in this experiment is obtained with a forcing period of 2.5 
days. Since the Doppler-shifted frequency may be written 
as to = k(c* - U Q ) = a) - kU o , this shorter period actually 
gives approximately the same Doppler-shifted frequency as 
that of Experiment I (where the period is 3.25 days with 
U Q = 0). Figure 25 shows the ip 1 field for day 70. 

Based on the results of the simple analytical model 
[Eq. (6.8)], we expect that the constant advective velocity 
U 0 will significantly reduce the local spatial growth rate 




Figure 26 shows the spatial growth rates from the parallel 
flow model (m) and from the complete numerical model (m*). 
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Experiment II. The field at t = 70 days. 
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Fig. 26. Experiment II. Spatial growth rates, m(x) from the parallel flow 
model, and m*(x) from the complete numerical model. 



Comparing this figure with Fig. 23, the reduction in both m 
and m* from those of Experiment I is clear. This reduction 
is evidently due to the advective velocity which causes the 
disturbances of Experiment II to move approximately 
-5 m sec" 1 faster than in Experiment I. Thus, the waves 
have less time to grow or damp in response to the local 
stability per unit distance traveled. The smaller spatial 
growth rate obviously affects the amplitude and spatial 
gradients of \p 1 . Figure 27 shows the wave packet envelope 
< ip 1 C x ) > of the complete numerical solution for Experiment II. 
Figure 28 shows this envelope along with the analytical model 
envelope F(x) using either x Q = -375 km or x Q = 3750 km as 
the initial point (no discernible difference is noted between 
the use of these initial points). Both envelopes are, as 
expected, significantly smaller than those of Experiment I. 

In addition, there is less variation in the wavelength 
(3750 km < x < 4875 km) compared to Experiment I. This is 
also a result of the faster phase speed which causes the waves 
to have less time to adjust to the local u(x,y). 

Although the faster phase speed of the waves in 
Experiment II reduces the magnitude of growth and decay rates 
and therefore brings the spatial growth rate curves of the 
numerical model (m*) closer to that of the parallel flow 
model (m), we notice a stronger asymmetry in the m* curve 
with respect to x=0 in Fig. 26 compared to Experiment I 
(see Fig. 23). In addition to the maximum m* occurring 
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slightly further downstream of the maximum of m(x=0) than in 
Experiment I, the cross-over points of the m* and m curves 
exhibit more asymmetry in Experiment II than I. The first 
cross-over point occurs approximately at longitude x = 8500 km 
and x = 4000 km for Experiments I and II, respecti ve ly , while 
the second cross-over point occurs at approximately 
x = -12000 km where the spatial growth rates vanish. As in 
Experiment I, we interpret this asymmetry as being due to the 
tilt structure of the waves lagging behind the expected value 
from the local stability condition. The fact that this lag 
effect is enhanced by the faster phase speed of the waves in 
Experiment II can also be seen in Fig. 29 where the north- 
south phase tilt 9 * ( x , y ) of the wave disturbances for both 
experiments are compared at several longitudes. The unstable 
tilt in Experiment II is weaker upstream and stronger 
downstream of approximately x = -2000 km reflecting the 
stronger lag effect of the faster moving waves. 

Figure 30 compares the phase tilt of the numerical model 
solution, 9*, and that of the most unstable wavelength of 
the parallel flow model, e , for two longitudes (x = 0 and 
x = 3750 km). Comparing this figure with Fig. 24 of 
Experiment I, it is obvious that the apparent augmentation 
of the wave tilt due to the zonal variation of u(x,y) is now 
significantly reduced. This behavior may at least be 
partially explained by the lag effect caused by the advective 
current U 0 , which generally reduces the growth rate in the 
vicinity of the jet maximum. 
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Phase Angle 6 (Degrees) 
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Comparison of phase tilt of the complete numerical model, 6*-, between 
Experiment I (dashed) and Experiment II (solid) for 1 ong i tudes ( a ) 
x = 13125 km, (b) x = 5625 km, (c) x - 0, (d) x = -5625 km, (e) 
x = -13125 km, (f) x = -19500 km, and (g) x = -24375 km. 
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Experiment II. Comparison of phase tilt of the complete 
numerical model, o*, and of the most unstable wavelength 
of the parallel flow model, o , for longitudes: (a) x = 
and (b) x = 3750 km. 



The lower portion of Fig. 31 shows the wavelength 
measured from the numerical solution at y = -750 km. In 
the upper portion of the diagram are the phase velocities at 
y = -750 km, c* from the complete numerical model and c from 
the parallel flow model based on the corresponding numerical 
solution wavelength. The two phase speeds agree very closely 
upstream of the jet maximum while c* is smaller than c down- 
stream of this longitude. Similarly, Fig. 32 contains the 
wavelength and phase speeds for latitude y=0. At this 
latitude, the two phase speeds have similar behavior through- 
out the length of domain with c* being approximately 
1 m sec - "' slower than c. Figure 33 shows the wavelength 
curves [ L ( x ) ] for latitudes y = 0 and y = -750 km. At both 
latitudes the maximum wavelength occurs near the jet maximum, 
but the wavelength (and therefore phase velocity) at 
y = -750 km is larger upstream and smaller downstream of 
approximately x = -5000 km. As in Experiment I, these 
differences in c*(x) and L(x) between the two latitudes are 
consistent with the observed tilt behavior of the waves. 

Other than the foregoing mentioned differences, the 

basic behavior of the waves in Experiment II is similar to 

that of Experiment I. This includes the tilt of the waves 

and its reversal from that of growth to that of damping 

downstream of <i>'> v » the occurrence of <<!>'> , v a consider- 

max max 

able distance downstream of the jet maximum, the general 
shape of the wave packet envelope, and the longitudinal 
variation of the energy balance (Fig. 34). 
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X ( 1 000 KM) + 

32. Experiment II. Latitude y = 0. (a) Phase velocities c(x) from 

the parallel flow model, and c*(x) from the complete numerical 
model. (b) Wavelength L(x). 
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33. Experiment II. Wavelength L(x) for latitudes y=0 and y = -750 km. 
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X ( 1 000 KM) + 

Experiment II. Disturbance energy balance: (a) Reynold's stress term. 
<( u ' v ' 3u/8y )>; (b) energy transport term <( -uak ' /3x)>; (c) dissipation 

term, <(-2D f k')>; and (d) pressure work term, <( - -•^-(u'p')>. 



VIII. SUMMARY AND CONCLUSIONS 



In this work we studied the behavior of waves moving 
into and out of a barotropi cal ly unstable mean wind field 
which varies both in x and y. This mean wind field roughly 
simulates the 200 mb easterly jet south of the Tibetan high 
during the Northern Hemisphere summer. A rectangular domain 
is used with a time periodic forcing on the inflow (east) 
boundary and a Sommerfeld radiation condition on the outflow 
(west) boundary. This allows the simulation of the propaga- 
tion of small amplitude waves through the easterly jet region. 
The vorticity equation on a beta plane is solved with the use 
of finite differences, and, when the boundary conditions are 
properly adjusted, the barotropic waves move smoothly across 
the region and out the western boundary. After a certain time 
period of numerical integration, the solution becomes 
periodic everywhere with the forcing frequency which is 
specified on the eastern boundary. As the waves move 
through the jet regime, they grow or decay spatially in 
response to the local stability of the mean flow, even 
though at each point the variation is purely harmonic. The 
numerical results are compared with a simple mechanistic 
analytical model which is developed using the local growth 
rate and phase velocity computed from the parallel flow 
theory . 
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Two experiments are carried out in this investigation. 
The parameters used in these experiments are the same, 
except that in Experiment II a constant advective velocity 
of U Q = - 5 m sec”^ is added to the mean flow in order to 
gain a better understanding of the adjustment process of a 
faster moving wave train. 

The basic behavior of the waves in both experiments is 
similar. This includes the following: 

1) The maximum amplitude of the waves occurs a 
considerable distance downstream from the most unstable part 
of the jet near where the local growth rate of the parallel 
flow theory vanishes. 

2) The tilt structure of the waves and its reversal 
from that of growth to that of damping occurs downstream of 
the maximum of the wave packet envelope in response to the 
local stability of the mean flow. 

3) The wavelength (and therefore phase velocity) 
increases slowly from the inflow boundary to a maximum value 
near the jet maximum and then decreases slowly downstream 
reaching a minimum value near the outflow region. 

4) The Reynold's stress term < ( u ' v ' au/ 3y ) > is the 
dominant source for disturbance kinetic energy and the 
advective term <(-u3k'/3x)> plays a significant role in 
transporting a considerable part of this energy downstream 
from the source region. 

5) Large spatial damping occurs in the outflow region 
due to the presence of the continuous spectrum modes. 
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In Experiment I where the mean wind vanishes toward 
large y, the maximum spatial growth rate in the numerical 
model is significantly larger than the spatial growth rate 
of the parallel flow model. Apparently the downstream 
variation of the mean wind augments the phase tilt of the 
waves thus causing an increased growth rate as compared with 
waves in a parallel flow. In Experiment II where the waves 
move westward faster due to the presence of a constant 
-5m sec"^ advective zonal wind, the spatial growth rate 
computed by the numerical model is brought closer to that 
given by the parallel flow theory, but the asymmetry of the 
numerical growth rate curve with respect to the jet maximum 
is enhanced. The reduction of spatial growth and decay rates 
is mostly due to the faster phase speed which causes the 
waves to have less time to grow or damp per unit distance 
traveled. The enhanced asymmetry is apparently due to a lag 
effect wherein the adjustment of the waves is delayed 
spatially. The latter effect is also reflected in the north- 
south phase structure of the waves and in the smaller varia- 
tion of the wavelength as the waves propagate downstream. 

The augmentation of the phase tilt by the downstream 
variation is also significantly reduced which may at least 
be partially explained by the lag effect. 

Recently Pedlosky (1976) used a simple analytical model 
to study the effect of the downstream variation of a baro- 
clinic current. The mean flow in his model changes abruptly 
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from a weakly unstable regime with a constant growth rate 
to a stable regime downstream. His model includes nonlinear 
effects but it does not produce strong dynamic damping due 
to the continuous spectrum modes. This is perhaps because 
he did not include the lateral shear of the mean flow. Thus, 
Pedlosky found that his disturbances remain substantially 
undiminished when they propagate from a locally unstable 
region to a locally stable region. 

The behavior of the wave disturbances in our study 
resemble certain aspects of the waves observed by 
Krishnamurti (1971a, b) along the South Asia easterly jet 
during the Northern Hemisphere summer, and the waves that 
were simulated by Colton (1973) in a nonlinear model. The 
wave disturbances entering upstream of the easterly jet 
regime grow and increase their speed and wavelength as they 
approach the jet maximum. They then achieve their maximum 
amplitude a considerable distance downstream of the jet 
maximum and subsequently decay as they leave the strong jet 
regime. 

The simple mechanistic analytical model developed in 
this study can crudely approximate the adjustment of the 
moving disturbance to the local stability of a variable mean 
flow within the dynamically unstable region. It is possible 
that this simple approach could be used to study the linear 
instability of other variable mean flows to wave disturbances. 
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Although this study employs linear equations, it is 
possible to estimate qualitatively the nonlinear effect of 
the waves on the mean flow. Since the Reynold's stress term 
is the dominant term in the energetics, most of the mean 
kinetic energy is removed by the waves on the downwind side 
of the jet maximum between where the wave amplitude is 
maximum and where the phase tilt of the waves is strongest. 
Clearly these synoptic scale waves will affect the amplitude 
and phase of the quas i -s tati onary planetary waves which 
combine to form the easterly jet over South Asia. These 
effects can be studied by including nonlinear effects in 
the present model. It may also be interesting to extend the 
present model to a 3-dimensional model to study the effect 
of vertical shear and barocl i ni ci ty . 
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APPENDIX 



Subroutine POISDD 

Subroutine POISDD is a direct method solver of the 
finite difference approximation to Poisson's equation 
(Sweet, 1971 ) : 

v 2 u(x ,y ) = - f (x,y ) (1 ) 

for a<x<b 
c<y<d 

and with Dirichlet boundary conditions 
u ( x , c ) = g i (x) 

u(x,d) = g£(x) a<_x<_b (2) 

u ( a ,y ) = 9 3 (y) c<y<b 
u(b,y) = g 4 (y) 

We define a grid on the rectangle {(x,y): a <x<_b , 

c<_y<_d} by selecting two positive integers M and N, such that 

AX = and AY = ^ (3) 

and where N must be a power of two, and the points 

x. = a + (i-l)Ax i=l, 2, . . . ., M+l 

1 (4) 

Yj - c + (j-l)AY j = l, 2 N + l, 

where i and j refer to the x and y grid points, respecti vely . 

The Poisson equation that is solved by POISDD is given 

by Eq. (2.10): 

V 2 ||^- = - J(t,V 2 ip') - JU',v 2 ij;) - - D f vV 
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with boundary conditions Eqs. (2.11), (2.14) and (2.16). The 
finite difference form of these equations are given by Eqs. 
(4.10), (4.14) and (4.17). Here T . ^ . is the finite differ- 

' iJ 

I 

ence representation for the tendency, dip /9t, and is the 
dependent variable in our Poisson equation. For our 
Numerical Model, M=106 and N=32; t is the current time, 
where t = £At i = 1,2.... and At in the time step. 

Our system of equations meets the requirements of 
subroutine POISDD, except for the outflow boundary condi- 
tion, Eq . (2.16). 



9 / 3 ^ 1 \ _ 9 / djjj 1 \ 

3 1 3 1 ci v * ^ • 



r 9X v 9t 
0 



( 6 ) 



This equation is the Sommerfeld radiation condition (Pearson, 
1974). Clearly, this is not the usual form of Dirichlet 
boundary condition. Therefore, in order to use POISDD in 
our numerical model, this subroutine had to be modified 
incorporating the radiation outflow condition. This was 
done by writing Eq. (6) in finite difference formr 



T t T t-At 

'u ' 1 1 , j 

At 



T t - T ^ 

2 , j 1 1 , j 

AX 



(7) 



where the finite difference symbols have their usual meaning 
in numerical weather prediction. Subscript i=l refers to 
the outflow boundary, and subscript i=2 refers to the 
column of grid points next to this boundary. We observe 
that this scheme is backward in time with upwind space 
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differencing. Equation (7) is conveniently rewritten such 
that 



yt-At 



To t • 



T 1 , j ’ 1 + c' At/AX + 1 + Ax/ ( C At) 

r o r o 

Since all the interior points of the domain must satisfy the 
non -homogeneous Poisson equation, we, then, use this constraint 
to. couple the radiation boundary condition [Eq. (8)] with the 
interior points, as follows: 



- T i!i , T 2 !.i+1 - 2T 2 !i + T 2 !i- 1 

ax^ Ay 2 



(9) 






The term f. . is the finite difference approximation to the 

I )J 

right side or non -homogeneous part of Eq . (5) or Eq . (2.10). 

2 

We multiply Eq . (9) by -(Ay) , and we obtain 
T t 2T ^ T ^ 

. t ii£j. . Iiii . T -t * 2T t . T t 

s s s 2 , j + 1 2 , j 2 , j -1 






2 2 

where s = ax /Ay 



We now combine Eq . (8) and Eq . (10), and we get 

T t 



' T 2, j + l + 



2 + - - 



1 



s(l + 



AX 



•) 



C At 
0 -J 



2 , j 



t-At 



3 , j + 
s 



T t _ - . 2 1 , ,i 

’ ' 2 , j + 1 " f 2,j Cyy + sLl+c" At/AxJ 



( 10 ) 



( 11 ) 
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By incorporating Eq. (11) in subroutine POISDD, we 
have, thus, coupled the Sommerfeld radiation boundary condi- 
tion to the interior domain. Also, by Eq . (8), we now have 
a complete set of Dirichlet boundary conditions. The finite 
difference Eqs. (4.10), (4.14) and (4.17), thus, form a 
linear system of equations of dimension (M-l ) x (N-l) for 
the unknowns T . ^ : 2<_i<M, 2<_j<N. This system is solved in 

* J 

POISDD by the Buneman variant of the CORF (cyclic odd-even 
reduction and factorization) algorithm. Buzbee et_ aj_. ( 1 970 ) 
gtves a complete mathematical description of the algorithm. 
Sweet (1971) gives a brief description of the subroutine. 

A function of two variables, B. . , is used in the 

■ 5 J 

numerical model code to represent a two-dimensional array 

which provides, on input to POISDD, values of the function 

f. . as well as the specified boundary conditions. On output 
1 5 J 

from POISDD, it provides the values of the approximation 

For our numerical model, we approximate the outflow 
boundary condition on input to POISDD as follows: 

T t-At 



T* 
i . J 



B 



1 J 



1 , j 



1+C At/Ax 
r o 



( 12 ) 



On output from POISDD, we correct Eq . (12) by 



T t _ n * + 2 ,_j _ 

T l,j B l,j l+Ax/(c~ At) 



(13) 
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